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Abstract 

3-D BEM formulations of unsteady lifting flow problems are paradigms of free boundary problems. The free boundary enters in 
the BEM formalism as the necessary spatial support of the velocity discontinuity, generated at the trailing edges of the lifting 
bodies. In the context of a BEM, this velocity discontinuity, alternatively termed a shear layer, is simulated by a surface dipole 
distribution or its equivalent singular vortex sheet. Free vortex sheets are inherently unstable and amenable to two well-known 
instabilities: the Kelvin-Helmholtz instability and the roll-up of their free edges. As a result, they show chaotic behavior with 
the passage of time. On the other hand, experimental evidence on flows around lift producing devices shows that the wake 
vorticity is organized in specific, problem dependent rollup patterns. In this paper we present a filtering technique, which by 
introducing artificial viscosity to our problem, suppresses the smaller scale instabilities, leaving the larger scale organized 
vortices to determine the rollup pattern of the free wake vorticity. Flow simulations are presented for a steadily advancing wing, 
a biomimetic wing at two different Strouhal numbers, a naval propeller at two different advance coefficients and a novel 
flexible oscillating duct propulsion device. The effect of the filtering parameters to the shear layer wake rollup pattern as well as 
to the calculated unsteady forces and moments is shown and discussed. 
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Introduction 

Boundary element methods (BEM) have been long applied for the solution of flow problems around systems of 
streamlined lifting or non-lifting bodies. In the excellent book of Katz & Plotkin (2001) [10], almost all aspects and 
the state of the art of modern BEM formalisms, for fluid flow problems, are presented and discussed. BEM 
formulation of flow problems around lifting bodies is a paradigm of a Tree boundary problem'. The free boundary 
enters as the spatial support of the tangential velocity discontinuity, generated at the trailing edges of the wings. In 
the context of a BEM, this tangential velocity discontinuity is simulated by a dipole distribution or its equivalent 
singular vortex sheet. Typical remedies for the determination of this free boundary are: (a) the prescribed wake 
shape method or (equivalently) the frozen wake model by which the geometry of the free boundary is somehow a 
priori selected, (b) the wake relaxation method by which an initial supporting geometry is assumed and a corrector 
method based on Helmholtz vortex theorems is applied to obtain the final geometry, and (c) the time stepping 
algorithm by which the supporting geometry is derived by applying the free shear layer kinematic and dynamic 
boundary conditions in a step by step basis, Politis (2004) [18], Wu et al. (2006) [35]. Remedies (a) and (b) are 
applicable in either steady flow problems or in unsteady flow problems, where the unsteady flow field does not 
affect considerably the free shear layer position. Remedy (c) is more general and can be applied without limitations 
to the solution of any steady or unsteady flow problem. Unfortunately, application of remedy (c) requires coping 
with the inherent free singular vortex sheet instabilities. More specifically, it is well known that a singular vortex 
sheet is amenable to chaotic behavior due to a Kelvin-Helmholtz instability. Furthermore, the edge of the sheet 
induces infinite velocities to itself, fact which is related with the tip vortex rollup. Finally, a singular vortex sheet 
induces a non-physical tangential fluid velocity jump in its vicinity, not-existent in a viscous environment. An 
extensive discussion of the free vortex sheet instabilities, their chaotic behavior and corresponding simulation 
methods can be found in the review by Sarpkaya (1989) [26]. The same subject is discussed in the books by Saffman 
(1992) [25], Marchioro & Pulvirenti (1994)[12], Bajer & Moffatt (2004)[2] and Wu et al. (2006)[35]. Furthermore, a lot 
of experimental and theoretical/numerical works have been done on the subject, such as: Winckelmans & Leonard 


1 



www.daoe-journal.org 


Development and Applications of Oceanic Engineering (DAOE) Volume 5 2016 


(1993) [33], Ramsey (1996) [22], Riley & Lowson (1998) [23] in connection with delta wings. Wood & Grace (2000) [34], 
Krasny et al. (2002)[11] in connection with the vortex blob numerical simulation method, Dong et al. (2006)[9], 
Muijres & Lentink (2007)[14], Parker et al. (2007-2008)[16-17], Borazjani & Sotiropoulos (2008)[3], Wang et al. 
(2008) [32] and Bohl & Koochesfahani (2009) [4], to mention some of them. 

Correct prediction of the dynamic evolution of the free shear layers is of high importance in many physical 
problems. For example, in biomimetic flows, the bird or fish wing performs large oscillatory (flapping) motions in 
comparison to its parallel translation, producing a highly deformed shear layer in the vicinity of the trailing edge, 
which affects the wing flow field. Similarly, in steady or unsteady naval propeller operation, the free shear layers 
are the main mechanisms of blade-to-blade interaction and introduction of the history effects. For higher loadings, 
the shear layers are limited in small regions behind the biomimetic wings or the propeller disk. As a result, a strong 
amalgamation of the free shear layers occurs and the problem of the mutual interaction of neighboring shear layers 
is superimposed to that of their self-action. 

In Politis (2004) [18], we have formulate and solve the problem of flow around a propeller performing a general 
unsteady motion. A potential based BEM has been used to represent the propeller blades while a time stepping 
algorithm has been used to simulate free shear layer dynamics. During the period 2002-2008, this formulation, 
initially developed to treat propeller problems, has been expanded to include the simulation of any 3-D unsteady 
incompressible non-viscous flow problem around an arbitrary system of interacting non-lifting/lifting, 
rigid/flexible streamlined bodies. By flexible we mean pseudo-flexible; that is body deformation is determined by 
the user and not through coupling with the elastic behavior. The developed code has also been enriched with two 
new non-linear Kutta conditions of the pressure type, in addition to the initially used Morino type Kutta. The 
corresponding formulation is discussed in Politis (2009, 2011) [19-20] where comparisons with experimentally 
measured forces for biomimetic systems and propellers are also presented. 

In this paper, we concentrate on the methodology used for the numerical determination of the free shear layer 
dynamics and more specifically, to the filtering technique, which suppresses the smaller scale instabilities leaving 
the larger scale dominant vortices to specify the dynamic evolution of the free shear layer surfaces. The method 
works successfully for any flow case, including that of large unsteady motions of interacting lifting rigid/flexible 
bodies. After a presentation of the method, we give a number of applications as follows: (a) a burst starting steadily 
advancing 3-D wing, (b) a biomimetic 3-D wing operating at two different Strouhal numbers, (c) a burst starting 
steadily advancing and rotating naval propeller operating at two different advance coefficients, and (d) a novel 
FOD (Flexible Oscillating Duct) propulsion system. The last application shows how our method can be used as a 
tool in analyzing and/or design of novel hydrodynamic systems based on the new Electro-active Polymer (EAP) 
materials, which can actively deform in a fluid environment, simulating the action of muscles, Cohen ed. (2004) [5]. 
In all cases, we show the effect of the filtering parameters to the obtained shear layer 3-D patterns. It is shown that 
for a range of the filtering parameters the free shear layers are organized in concrete 3-D space patterns, 
characteristic to the problem considered. We thus arrive at the notion of the 'attracting wake configuration' of a 
problem. Furthermore, calculated unsteady forces and moments are almost independent of the filtering parameters, 
even in chaotic shear layer pictures, as far as the shear layer geometry lies within the attracting configuration. 
Finally, it is shown that in certain flow problems, the old assumption of a frozen wake model results in unsteady 
forces and moments comparable to that using the more elaborate free wake model. 

Formulation 

We start with a brief introduction of the main theoretical & numerical aspects of the BEM, necessary for the 
description of the filtering technique. A more detailed discussion can be found in Politis (2011) [20]. 

Geometric Considerations 

We build the geometry of complex systems of bodies using surface patches. Each patch consists of a number of 
bilinear quadrilateral elements. Two types of patches are allowed: (i) lifting patches and (ii) non-lifting patches. By 
combining patches we build lifting and/or non-lifting bodies. For the former case, the user has to determine the line 
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of flow separation in the surface of each lifting patch. The total system surface at time t is denoted by SB(t ) . Total 
system Kutta strip surface at time t is denoted by SK(t ) . Total system free shear layer surface at time t, excluding 
Kutta strips, is denoted by SF(t), Politis (2011)[20]. 


Velocity and Potential Representation Theorems 

We use an inertial (build in earth surface) frame of reference for the definition of velocities (body or fluid). A 
corresponding coordinate system (assumed Cartesian-orthogonal) is denoted by OXYZ. 

As a result of the (known) unsteady motion of our system of bodies, in the region outer to SB(t)u SK(t)u SF(t ) , a 
velocity (perturbation) potential (/) exists which, at each time step, is expressible through its traces </>, V $ on the 
boundary points Q e SB(t ) u SK(t ) u SF(t ) . 

Introduce: 


1 r h-Vd> 1 r ,n-r 1ry 1 r n-r 1 r 

F(P) = —dS H 6——dS H u — -dS-\ 

J y \n ^ r 4 tc ^ r 4/r ^ 

SB(t ) 55(0 ' SK(t) ' ^ SF(t) 


n ' T 
JU — 


(1) 


H(P) = — f (n-V(zS)4^+— f (nxV^)x-CrfS+— f fx-^dS+ — f f 
Att r 4 jtt ^ r 4^r ^ r J 4^r ^ 

55(0 56(0 SK(t) ' 55(0 


' x — — dS — 


4tt 




dl xf 


( 2 ) 


L'(t) 


where: P is the evaluation point (or control point) for either F or H , h is a unit normal vector at the boundary 


QP 


, ju is the dipole 


integration point Q g SB(t)u SK(t)u SF(t) showing inside the flow region, r = QP , r = | 
intensity with support on SK(t)u SF(t) and y the corresponding (to ju) surface vorticity intensity given by: 

ju = f-(j)- (3) 

y — hxV p (4) 


Finally, L\t) (integration region of the last line integral in the right hand side of equation (2)) is the free part of the 
line bounding the free shear layers, defined by: L\t) = d(SK(t) u SF(t )) - ( SK(t ) n SPb)) . 

With the aid of relations(l), (2) representation theorems for (j), V(/> become: 
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V^P) = H{P) 
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>P e 5 , P + ” (?) 
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(6) 


(7) 


and similarly for SX + . In relation (7), the superscripts ( + “) denote the two sides of the free shear layer surfaces, 

while the unit normal n is directed from (-) to (+). Notice that some of the surface integrals for F or H in relations 
(6) or (7) contain strong surface singularities of the Cauchy type. Thus, their meaning is realizable only in the 
principal value sense, Mikhlin (1965) [13]. 


The Integral Equation 

Let v A denotes the (known) velocity of the boundary point A e SB(t) and h a unit vector normal to body surface at 
A with direction pointing into the flow region. Then, the no-entrance condition at A has the form: 
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V (j)-n = v A -h 

Substituting ( 8 ) to the first of equations ( 6 ) and using (1) we get: 



This is a second kind Fredholm type Cauchy singular boundary integral equation for the determination of ^ and 
ju on points of SB(t ) and SK(t ) respectively. In the right hand side of (9), the first term is a known integral (as far 
as the motion of the system of bodies is known). The second term is an integral over the free boundary SF(t) 
known from the solution of the problem at previous time steps. The unknowns in the left hand side of (9) are the 
potential (f> on SB it ) and the dipole intensity ju on SK(t) . For their determination, the additional required 
condition is the Kutta condition at the separation lines (trailing edge in case of a wing flow w/o separation). 


The Kutta Condition 


Let the point A e SB it ) . Let — 

dt 


denote the time derivative for an observer built to the point A of the moving body 


and let v A denote the known velocity of A . Then unsteady Bernoulli equation takes the form: 


P-P x 


dcf 

dt 


4( v ^) 2 +k 2 


( 10 ) 


According to a pressure type Kutta condition, as we approach the trailing edge point from either pressure side 
(superscript +) or suction side (superscript -), the pressure should be continuous, i.e.: p + = p ~ . Using(lO), this 
becomes a quadratic (nonlinear) relation between ^ + , V^' . Assuming steady linearized flow, Bernoulli 
equation degenerates to the famous Morino condition: jtt amoachingL(x) = (^ + -(/> ") approaching l© which is a linear equation in 

from shear layer from body points 

(j), jLi . Either form of Kutta conditions can be handled by our code. 


Shear layer dynamics. Kinematic and dynamic conditions on a free vortex sheet expressed in terms of the dipole 
intensity of the sheet result in the following equation, Politis (2004)[18], Wu et.al.(2006)[35]: 



(ii) 


where D / Dt denotes a material derivative for ju based on the mean perturbation velocity > v < of the shear 
layer: 


V dP + V ^ 

> v <= — — = H(P) 


( 12 ) 


or in developed form: 


> v(P) <= 




r 1 r r 1 r 

x— ds +— f fx— ds-— f JU 

r ^sn» r 4 k Ao 


dl xr 


(13) 


Equation (11) informs us that the dipole surface SF(t ) , supporting the dipoles with intensity is travelling 

with velocity > v < , where <^,77 denotes a set of curvilinear surface coordinates for the points on SF{t ) . Thus, if a 
ju surface exists at time t , we can find its new position, at time t + dt , by deforming it by > v < • dt . 


Calculation of Forces, Moments and Power 

Let Fit) denote the total force exerted by the fluid to the body at A e SB(t ) . This force is a sum of pressure forces 
normal to the body (denoted by P(t) ) and viscous forces tangential to body (denoted by D(t) ). Thus: 
F{t) = Pit) + Dit ) . Pressure forces can be calculated from equation (10) while drag forces can be found assuming a 
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local at A viscous drag coefficient. Assuming a reference point for moments r M (t) , the total instantaneous system 
force/moment can be calculated by the formulas: 

J J r M (t)x F(t)dS (14) 

SB(t) SB(t ) 

The last formulas can be applied to the whole system (use SB(t ) as integration region), or to one or more bodies of 
our configuration (use a subset of SB(t) as integration region). In propulsion problems, there is always a preferable 
instantaneous direction in which the system moves. Let V(Othe instantaneous velocity of the system along this 
direction. Then the instantaneous useful power of the system is given by: EHP(t) = V(t)- F m (t) ( EHP : Effective 
Horse Power). The net instantaneous power interchange of our system with the environment is given by: 
J F(t) • v A {t) • dS . Then the power provided to the system, termed DHP (Delivered Horse Power), is given by: 

SB(t) 

DHP{t)= J F(t) ■ v A (t) ■ dS - EHP(t) (15) 

SB {t) 


The ratio EHP(t) / DHP(t) defines the instantaneous efficiency of our system. Notice that the previous formulas can 
be applied equally well to either rigid or flexible bodies. 


Discretization and Solution 

Subdivide SB(t) into Af 5 elements. Subdivide SK(t) into A/^ elements. Subdivide SF(f)into A/^ elements. Four node 
quadrilateral elements have been used for the subdivision of body and shear layer boundaries. Assume piecewise 
constant (j) and n-V(/> for all elements on SB (t) . Assume piecewise constant ju for all elements on SK(t) u SF(t ) . 

Denote these constant values by $,cr.(= (n- V (/>).), ju t where the range of index (i) is adapted accordingly. With the 
aid of the previous assumptions/notation integral, equation (9) becomes: 

Z B ‘Ji~ Z Z A .&+ Z B uBj ( 16 ) 

Z i=i,N B j=\,N K j=i,N B j=l,N P 


J_r JS_ B 1 t n(Q)-QP,dS 

^\\QPi r u 4 *l \qp\ i 


(17) 


and Ej denotes the surface of the j th element from either SB(f),SK(t\SF(t) and P denotes the i th control point 
(centroid of E t ) on SB(t) . Relation (16) applied at the N B centroids of the body elements, gives N B linear 
equations for the determination of the (j) j , // ; . The N K additional equations required for the calculation of /Uj are 
taken from the satisfaction of the Kutta condition on SK(t)nSB(t) (i.e. trailing edges). In our code we implement 
three alternative forms of the Kutta condition: (i) The first alternative consists of a linear Morino condition, which 
in discretized form becomes: = </>. + - (/>._ ( i + and i denote element numbers on body, neighboring to trailing 

edge from different sides and i denotes element number on kutta strip, neighboring to the same point of the 
trailing edge), (ii) The second alternative consists of a nonlinear pressure type Kutta in the form of equation 
p + = p~ with the additional conditions that the downstream length of the Kutta-Strip element is calculated 
according to the total velocity there and the kutta-strip element is either tangential to blade face or blade back 
according to the local value of circulation at this span- wise position, (iii) The third alternative consists of a 'mixed 
type Kutta' i.e. partly Morino and partly Pressure type. We have found that at higher loadings, a pressure type 
Kutta at blade tips is too strong and occasionally leads to a destruction of the shear layer geometry at those regions. 
Thus, we have decided to introduce a third alternative, the "mixed Kutta", i.e. a Morino condition at the tips of the 
separation line and a pressure type Kutta at all other points. If alternatives (ii) or (iii) have been selected, the 
resulting system of equations is nonlinear and it is solved by using a Newton iteration method, with starting value 
taken from a Morino type Kutta (first alternative). 
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With the previous discretization scheme in mind, the discretized formula for the calculation of > v < at points on 
SF(t) , needed for the application of the shear layer deformation scheme, becomes: 

>v,.<= I C,.< 7 . + X D,J ]+ X D ijMj + X D iiMj (18) 

j=l,N B j=l,N k j=l,N F 

J HQ) x QPi dl ^ j a un - t tangent vector to dE j ) (19) 

d Ej \QP t \ 


where: 


j=hN B 


' if 


QP,dS 


\qp\ 


" , 


For the derivation of formula for D j j the equivalency of a constant dipole element with a line vortex at its 
boundary has been used, Politis (2004) [18]. Furthermore, for a bilinear boundary element EV, its boundary dE . 

consists of four straight lines. Thus, the corresponding line integral in equation (19) can be calculated analytically. 
In detail: let dEj = A5u BC u CDkj DA ( A,E,C,E> denote the nodes of element E j ). Then: 

D. . = D ; . + D . + D. . + /3 . . The contribution D. . to D . from the side AB (and similarly from sides: 

l ,J hj \ab l ’J\BC l ’ J \ CD l ’ J \ DE l,J \AB l,J V J 

BC , CD , DA ), to the control point P , is given by the relation: 


D 


1 

4zr/z 


(coscir + COS /?) 


PAxPB 

PAxPB 


where a = j£AB,AP , j3 = £BA,BP and /z denotes the normal distance from control point P to the line AP . 


(20) 


Solution of our problem is implemented by a time stepping algorithm as follows: 

At each time step: 

a. Find the next geometric position of the system of bodies. 

b. Generate corresponding Kutta strips, for the case of lifting bodies, introducing thus the extra unknowns 
required for the Kutta condition satisfaction. 

c. Solve the system consisting of the "no-entrance", equation(16), and "Kutta" conditions. In case of pressure 
type Kutta, a Newton iteration is used at this step. 

d. Deform the free shear layers to their new positions, by applying a special filtering technique (to be 
discussed later) to calculate > v < , equation(18). 

e. Output results (pressures, forces, velocities, position of free shear layers) at this time step, in a form of 
TECPLOT compatible files. 

f. Proceed to the next time step and repeat the calculation starting from step (a). 


The Filtering Technique 

From the phenomenological point of view, the physical mechanisms of either the self-action of a free shear layer or 
the mutual interaction between neighboring shear layers, are mainly viscous and as such are excluded from the 
modeling capabilities of a rigorous BEM formulation based on the representation theorems for potential and 
velocity. 

More specifically, determination of either the shear layer self-motion or its interaction with the neighboring shear 
layers requires evaluation of the velocity > v < , equation(18), at control points P e SF(t ) . In our formulation we use 
as control points, the grid nodes of the moving grid representing the free shear layer at time t . Considering the 
terms in the right hand side of formula(18), we observe that only the last two contribute to > v < in a singular 
manner. This singular manner is explicitly expressed by relation (20). A heuristic proposal for the elimination of 
the singularity is offered by either the / Lamb-Oseen / or the 'Burgers' vortices, which are special exact solutions of 
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the Navier-Stokes equations. In both cases, the azimuthal velocity U e of a vortex filament is given by the equation, 
Wu et al (2006) [35]: 


r - c (-T 

u e =—{ l-e {r) ),p = 2 

2 nr 


(21) 


where: T denotes the intensity of the vortex filament, r the radial distance of the control point from the core of the 

- c (-T i — 

filament, c is a constant characterizing the range of action of the term ( l-e yR ) and R- 2-Jv-t is a time 
dependent characteristic length used to make the action distance r non-dimensional ( V denotes the kinematic 
viscosity of the fluid and t the age of the vortex). 

Comparing formula (21) with formula (20) applied to the case of an infinite singular vortex line 
(a = 0° = 0° => cos a + cos j3 = 2 ), we observe that they are similar except for the multiplicative factor: 


m(c, p,r / R) = ( l-e ^ ), p = 2 


( 22 ) 


which 'mollifies' the singularity at r = 0 , substituting it by a viscous vortex core of finite diameter, controlled by 

the parameters: c,p,R. Rearranging the terms in (22) and taking the logarithms we get: c = -ln(l-ra)/ 

Assuming a characteristic intensity m = 0.5 atr/R = l, the mollifier constant C can readily been calculated as: 
c = -ln(0.5) = 0.69314718 independent of p . In this case, the mollifier function is uniquely defined by its radius R 
and its power p ( p = 2 for vortex filaments). 

Formula(21), for r = In and preselected^ = 0.69314718, p = 2 ,R = 1, is shown in figure 1. Formula (22) is shown in 
figure 2 for three different values of the parameter p = 2,3,4 and withe = 0.69314718 . 



Notice that mollifier functions have already been introduced in vortex particle flow simulation methods, 
Winckelmans & Leonard (1993) [33], Krasny et. al. (2002) [11]. 



FIGURE 1. AZIMUTHAL VELOCITY OF A LAMB-OSEEN VORTEX. T = 2tt,c = 0.69314718, p = 2,R = l 
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r/R 


FIGURE 2. MOLLIFIER FUNCTION, c = 0.69314718, p = 2,3,4 
With the aid of the mollifier function (22), formula (20) becomes: 

: -c{—\ 1 PA x PP 

D. . =(l-e [Rj ) (cos cc + cos fF) , „ __, ,/? = 2 (23) 

Ub 4tt/z PAxP# 

where ^ denotes the mollified equivalent to D j j . Considering formula (18) we observe that the D Uj is met at 
the second, third and fourth terms in its right hand side. Since the control point for >V,.< belongs to SF{t ) , only the 

third and fourth terms result in singular behavior and thus have to be mollified (i.e. use D iJ | instead of D Uj ). This 

observation is correct for all control points on SF{t) with the exception of the region of SF(t) in the vicinity of the 
Kutta strip (i.e. adjacent to the trailing edges). In this region, the induced velocity is strongly affected by both the 
free shear layer and the body elements. By applying a filtering process solely to the free shear layer dipole elements, 
an 'asymmetry' in the induced velocities in the vicinity of the Kutta strip can occur, which can trigger the 
breakdown of the shear layer geometry at those points. This danger is more pronounced at the first few time steps 
of a burst starting simulation, where a strong starting vortex is created in the vicinity of the trailing edges. To avoid 
the numerical instability, the second term in the right hand side of equation (18) (i.e. the induced velocity 
contribution from body potential) has to be mollified as well. Furthermore, for the same symmetry reasons, the 
mollifier function has to be smooth (or at least continues) in its extended support: SB(t)u SK(t)u SF(t ) . In our 
formulation, we satisfy this condition by taking the body mollifier radius R to be constant for all body elements 
with a value equal to the shear layer mollifier radius R at the trailing edge region. 

It is a well-known fact that viscous diffusion increases the diameter of the vortex core of a viscous vortex filament. 
This is the physical content of the relation R = 2-Jv-t introduced previously. Thus, older vortices have larger 
viscous diameters (or action ranges). According to the time stepping algorithm, a new span- wise strip of ring 
vortices is added to each shear layer at each time step. Thus, different span- wise vortex strips have different ages. 
In our formulation, we take heuristically this effect into consideration, by introducing age dependent mollifier radii 
R . This plays a crucial role in delaying the onset of chaos in the older parts of the free shear layers, especially in 
cases with higher hydrodynamic loadings. A measure of the age of a span- wise strip can be provided by the 
number of strips in-between this strip and the trailing edge. From the experience gained from the practical 
applications of the method, we know that proper values for R are: 0.5% to 15% of the diameter of the convex hull 
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(denoted by D CH ), containing our system of bodies, depending from the hydrodynamic loading, the body 
geometric details and the motion details of our system. 

In the case of a close contact of the freely moving shear layer with a body surface (for example in a kaplan type 
ducted propeller, the blade shear layers pass very close to the duct surface) numerical instabilities may occur due 
to the zeroing of the denominator in the respective formulas for the calculation of velocities and/or potentials, 
equations (18) and (16), which can lead to an intrusion of the shear layer inside the body. To avoid such an instance, 
the developed code contains an algorithm by which user controlled 'elastic shields' are applied to the bodies, 
avoiding thus the intrusion of a shear layer node inside the body. This shield methodology is planned to be 
discussed in a future paper. 

Applications 

The Shear Layer Geometry of a Burst (Step) Starting Wing 

The first application concerns the flow around a burst starting wing. The wing has a span s = 2m , a chord at mid- 
span c = hn and a chord at tip c tip = 0.3m. The wing has an X-skewback of 30 degrees, measured at the c 13 span- 

wise line from the leading edge. A NACA 0012 section has been used. Calculations have been performed for an 
advance velocity V = 5m/ sat an angle of attack a = 20°, using a blade grid with 40 span- wise elements and 16 
chord-wise elements for either face or back. We use an equally spaced span-wise and a cosine spaced chord-wise 
subdivision. Total simulation time is 1.0s subdivided in 84 time steps, with each time step equal to St = 0.0119s . 
The previous selections result in an X-extension of the free shear layer, approximately equal to 5 chord lengths. 
Furthermore, D CH = 2m for this example. A mixed type Kutta condition has been used in this simulation and all 
other simulations except the last (FOD), in which a pressure type Kutta has been used. 

Simulations have been performed for the following wake models (WM): (a) a free WM with mollifier radius 
R / D ch = 0.005 , (b) a free WM with mollifier radius R / D CH = 0.07 , (c) a free WM using a time dependent mollifier 
radius R(t) , and (d) a frozen WM. In case (c) the mollifier radius varies linearly between time steps, with the 
following break points ( its denotes the order of the shear layer span-wise strip with regard to the trailing edge): 
(i its,R(t )/ D ch ) = (1,0.005), (21, 0.03), (84, 0.1) . In case (d), the support of the free shear layer dipoles coincides with the 
trace in time of the trailing edge (i.e. located at the undisturbed flow position). 

Figure 3 shows the free shear layer geometry for WM (a) colored with the dipole intensity//. Since y - hxV/u , 
equation (4), the constant- // lines, shown in the figure (and all similar figures to follow), are tangent to the surface 
vorticity intensity vector y . Thus, figure 3 represents a characteristic quadrilateral ring vortex with three free sides, 
consisting of the two tip vortices and the starting vortex, and a bound side coinciding with the wing body. As we 
move downstream from the wing trailing edge, the tip vortex roll-up gradually develops, while approaching the 
junction point with the starting vortex instabilities appear. By increasing the mollifier radius to R/ D CH = 0.07 , 
WM(b), the tip vortex instability is suppressed but with a similar effect to the tip vortex roll-up, figure 4. To allow 
a simultaneous development of the tip vortex roll-up with a suppression of the starting vortex instability, a time 
dependent filter can be used. Corresponding results are shown in figure 5. With the time dependent filter, WM(c), 
the starting vortex instability has successfully been limited, without affecting the tip vortex rollup at the region 
near the wing trailing edges. 

Comparing the 3-D views shown on figures 3, 4 and 5, we conclude that irrespective of the selected filters, a 
definite shear layer 3-D pattern exists in the form of a 3-D ring vortex. A better sense of the effect of the filters to 
the shear layer pattern can be obtained, if we compare the slices of the shear layers with a plane normal to the X 
axis for different mollifier radii. Figure 6 shows such a comparison for the plane normal to the X-axis at the point 
x = —2m (with regard to the global inertia coordinate system). Solid line refers to WM (a), dashed line refers to 
WM (b) and dashed-dot line to WM(c). 

Attempts to simulate the shear layer geometry of a steadily advancing wing, using singular vortex sheets, has been 
made in the past by other researchers with results similar to that presented here. We refer to the works of Ramsey 
(1996) [22] and Wood & Grace (2000)[34]. 
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Figure 7 presents comparisons of the vertical (along Y-axis) force, divided by fluid density p , obtained using either 
of WM (a), (b), (c) or (d). The curves in figure 7 are labeled as follows: WM (a): 0.005, WM (b): 0.07, WM(c): Rif ) , 
WM(d): frozen . The results for the vertical force are practically the same for WM's: (a), (b) and (c). Those 
predictions differ from the corresponding frozen wake model by about 0.6% . Thus, a frozen or a generalized 
wake model can be an attractive alternative to the more elaborate free wake model, for predicting forces in a 
steadily advancing wing. This explains the success of the traditional lifting line or lifting surface theories in 
predicting the steady lift for wings using either a frozen or a generalized wake model. 



FIGURE 3. FREE SHEAR LAYER OF A BURST STARTING WING FOR WM (A). WING COLORED WITH POTENTIAL INTENSITY. SHEAR 

LAYER COLORED WITH DIPOLE INTENSITY. 



FIGURE 4. FREE SHEAR LAYER OF A BURST STARTING WING FOR WM (B). 

TIP VORTEX ROLL-UP PARTIALLY SUPPRESSED. STARTING VORTEX INSTABILITY SUPPRESSED. 
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FIGURE 5. FREE SHEAR LAYER OF A B URST STARTING WING FOR WM (C). 
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FIGURE 6. SLICE OF THE FREE SHEAR LAYER AT x = -2m .SOLID LINE: WM(A), DASHED LINE: WM(B), DASHED-DOT LINE: WM(C) 



FIGURE 7. VERTICAL FORCE OF A BURST STARTING WING. PREDICTIONS WITH DIFFERENT WM'S 
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The Shear Layer Geometry of a Biomimetic Wing 

Biomimetic wings are the subject of extensive investigations, since they are ideally suited in converting 
environmental (atmospheric or sea wave) flow energy to useful thrust, succeeding thus propulsive efficiencies well 
over one. An extensive collection of research works in the area can be found in the book of Taylor et al. (Eds.) (2010) 
[30]. The subject of aerodynamics of natural flyers is discussed in a recent book by Shyy et al. (2008) [27]. Reviews 
on the subject can be found in the papers of Triantafyllou et al. (2000) [31] and Rozhdestvensky & Ryzhov (2003) 
[24]. Experimental and numerical simulations of various biomimetic flow problems can be found in: Dong et al. 
(2006) [9], Muijres & Lentink (2007)[14], Parker et al. (2007-2008)[16-17], Wang et al. (2008)[32], Borazjani & 
Sotiropoulos (2008)[3], Bohl & Koochesfahani (2009)[4], Andro & Jacquin (2009)[1], to mention some of the most 
recent works on the subject. 


Due to the obvious importance of biomimetic flow problems, we have selected it as our next application in a 
biomimetic wing. In the study of natural flyers and swimmers in cruising condition, it is found that the Strouhal 
number (which can be considered as a measure of their hydrodynamic loading), defined by: 


St = 


z-f-K 

V 


(24) 


is within the range 0.2 < St < 0.4 , Triantafyllou et al. (2000)[31], Rozhdestvensky & Ryzhov (2003)[24], Taylor et al. 
(2003) [29]. In relation (24), / denotes the common frequency of heaving or pitching, /^denotes a characteristic 
heave amplitude and V the advance velocity. On the other hand, man-made biomimetic wings can be designed to 
operate in much higher hydrodynamic loadings (or Strouhal numbers). For those reasons, we have decided to 
present flow simulations at two different Strouhal numbers: St = 0.5 and St = 1.0 . In both simulations, all data have 
been kept constant except frequency / which has been doubled for the second case. 


Case with St = 0.5 

The selected wing geometry and BEM discretization is the same as that used in the previous example of a burst 
starting wing. The biomimetic wing advances with a velocity V = 5.0 mis, performs a heaving oscillation with 
amplitude =0.5 m, a pitching oscillation with amplitude 25.0 deg and a phase angle of 90.0 deg (with respect to 
heaving) and a common (for both heave and pitch oscillations) frequency / = 2.5 rps . The pitching axis has been 
selected at 33% from the leading edge at mid-span. The total simulation time is 0.85* or two periods. 60 time steps 
have been used per period or 120 time steps total. The previous data resulted in a Strouhal number: St = 0.5 . 

Figure 8 shows the free shear layer geometry for a mollifier radius R/ D CH =0.01 (where D CH =2 m as previously), 
colored with the surface dipole intensity ju . An extensive instability of the shear layer is observed: (a) at its edges, 
and (b) at its span-wise strips which correspond to a sign change of the sectional angle of attack (defined with 
respect to the total undisturbed fluid velocity). By increasing the mollifier radius to R! D CH =0.07, the instability is 
suppressed in the major part of the free shear layer in figure 9. From figures 8 and 9, we observe that irrespective of 
the filter used, the shear layers occupy a certain region of the 3-D space with a well-recognized pattern. For a better 
confirmation of this conclusion, we compare the slices of the shear layers with planes normal to either the X or the 
Z axis in figures 10 and 11. Deformation of the shear layers can be used as a visualizer of the large scale vortices 
created by the motion of a biomimetic wing. Using this deformation, we conclude that the wake of a biomimetic 
wing consists of a train of oblique vortex rings which produce a repeated S-shaped 3-D wake pattern (reverse 
Karman Vortex Street). We can arrive at the same conclusion if we use the constant dipole intensity lines (used to 
color the shear layers) in figures 8 and 9, to visualize the surface vortex lines. This S-shaped attracting 
configuration in figure 10 is characteristic for all thrust producing biomimetic wings and has been experimentally 
observed by many authors, Triantafyllou et al. (2000)[31], Muijres & Lentink (2007)[14], Parker et al. (2007-2008) [16- 
17], Zhang & Zheng (2009)[36]. 

Figures 12 and 13 present comparisons of the unsteady biomimetic wing forces obtained using a free WM with that 
obtained with a frozen WM. For the free WM, the following mollifier radii have been used: R/D CH =0.01 , 
R / D ch = 0.07 . Curve labeling follows that introduced in the example of a burst starting wing. It is shown that the 
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free WM results in practically the same X- and Y- forces, irrespective of the selected R / D CH . Those forces differ 
from that obtained using a frozen WM only slightly. For example, at t = 0.62s, which is near the fx / p, fy / p 
maximum values, the difference between the frozen wake and the free wake thrust force fx I p is less than 1 % . The 
corresponding difference for fy I p is 2-3% . Notice finally that according to our sign conventions, the thrust force 
in figure 12 is negative. This is a result of the motion of the wing which has been selected in the negative direction 
(with respect to the inertia coordinate system XYZ and to the fact that we have define positive 'forces' as that 
exerted to the body by the surrounding fluid). 



FIGURE 8. FREE SHEAR LAYER OF A BIOMIMETIC WING AT t = 0.8s WITH EXTENSIVE VORTEX INSTABILITIES. WING COLORED 
WITH POTENTIAL INTENSITY. SHEAR LAYER COLORED WITH DIPOLE INTENSITY. St = 0.5, R / D CH = 0.01 



FIGURE 9. FREE SHEAR LAYER OF A BIOMIMETIC WING AT t = 0.8s WITH INSTABILITIES SUPPRESSED IN THE MAJOR PART OF THE 

FREE SHEAR LAYER. St = 0.5, R/D CH = 0.07 
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FIGURE 10. SLICE OF THE FREE SHEAR LAYER AT 
t = 0.85, z = 0.0m . REVERSE KARMAN VORTEX STREET WAKE 
PATTERN. HAIR LINE: R/D CH = 0.01 , BOLD LINE: 
R/D ch =0.01 



FIGURE 11. SLICE OF THE FREE SHEAR LAYER AT 
t = 0.8 5, x = -1.8m .HAIR LINE: R/D CH = 0.01 , BOLD LINE: 
R/D ch =0.07 



FIGURE 12. THRUST (X-AXIS) FORCE OF THE BIOMIMETIC FIGURE 13. VERTICAL (Y-AXIS) FORCE OF THE BIOMIMETIC 

WING - PREDICTIONS WITH DIFFERENT WAKE MODELS. WING - PREDICTIONS WITH DIFFERENT WAKE MODELS 


Case with St = 1.0 

This example is similar to the previous one, except frequency, which has been doubled i.e. / = 5.0 rps . The total 
simulation time for this case is 0.4s or two periods. 60 time steps have been used per period or 120 time steps total. 

Figure 14 shows the free shear layer geometry for a mollifier radius R/ D CH - 0.01 . We observe that the instabilities 
shown for the lower Strouhal in figure 8 have been strengthened. By changing the mollifier radius to Rf D CH = 0.07, 
the instability is again eliminated in the major part of the free shear layer in figure 15. From figures 14 and 15, we 
(again) observe that, irrespective of the filter used, the shear layers show a well-recognized pattern. This pattern is 
also shown in the slices of the shear layers with the planes normal to either the X or the Z axis in figures 16 and 17. 
Using the deformation of the shear layers to visualize the large scale vortices, we arrive (again) to the repeated-S 
shaped wake pattern or reverse Karman Vortex Street. Comparing figure 9 with figure 15, we observe that by 
increasing the Strouhal number (and consequently the hydrodynamic loading) the train of the oblique vortex rings 
has been limited to a region closer to the wing body. 
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Figures 18 and 19 present comparisons of the unsteady biomimetic wing forces obtained using a free WM with 
either R / D CH = 0.01 or R/ D CH = 0.07 , with that obtained with a frozen WM. It is shown that the free WM results in 
practically the same X- and Y- forces, irrespective of the selected R / D CH . Those forces differ from that obtained 
using a frozen WM. For example, at t = 0.32s, which is near the fa / p, fy / p maximum values, the difference 
between the frozen wake and the free wake thrust force fa/ pis 6-1% . The corresponding difference for fy/ pis 
19 — 20% . Thus, by increasing the Strouhal number (equivalently the hydrodynamic loading), the differences 
between the frozen wake and the free wake become more pronounced. This is an expected consequence of the 
shear layer limitation in a region closer to the body surface as the Strouhal increases. Since birds and fishes operate 
at: 0.2 < St < 0.4 , it seems that the simpler frozen WM can be used in force predictions instead of the more elaborate 
free WM. 



FIGURE 14. FREE SHEAR LAYER OF A BIOMIMETIC WING AT t = 0.4s WITH EXTENSIVE VORTEX INSTABILITIES. St = 1.0, R / D CH = 0.01 



FIGURE 15. FREE SHEAR LAYER OF A BIOMIMETIC WING AT t = 0.4 s WITH INSTABILITIES SUPPRESSED IN THE MAJOR PART OF 

THE FREE SHEAR LAYER. St = 1.0, R/ D CH = 0.07 
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FIGURE 16. SLICE OF THE FREE SHEAR LAYER AT 
t = 0.45, z = 0.0m . REVERSE KARMAN VORTEX STREET WAKE 
PATTERN. HAIR LINE: R/D CH = 0.01 , BOLD LINE: 
R/D ch =0.07. 



time(s) 


FIGURE 18. THRUST (X-AXIS) FORCE OF THE BIOMIMETIC 
WING - PREDICTIONS WITH DIFFERENT WAKE MODELS. 




time(s) 


FIGURE 19. VERTICAL (Y-AXIS) FORCE OF THE BIOMIMETIC 
WING - PREDICTIONS WITH DIFFERENT WAKE MODELS. 


The Shear Layer Geometry of a Burst Starting Naval Propeller 

Our next application concerns estimation of the steady performance of burst starting naval propeller. The propeller 
is 3-bladed with a diameter D = 5.1m , a pitch ratio P/D = 1.1, a blade area ratio A E / \ = 0.55 and its detailed 
characteristics are similar to those of a B-series propeller, Oosterveld & Oossannen (1975) [15], but with slightly 
increased tip-skew. Calculations have been performed using a blade grid of 24 radial elements and 16 chord-wise 
elements for either face or back. Propeller hub is included in our BEM modeling. 

Calculations are performed at two propeller advance coefficients: j = _f_ = 02 and/ = 0.7 . The first concerns 

f-D 

operation at a heavily loaded point, while the second at a point with a lighter load, closer to propeller design point. 
Case with J = 0.2 

The propeller moves with a velocity of advance V = 2.28m /s and rotates with a frequency / = 2 rps or 120 rpm . Flow 
simulation lasts for two full propeller revolutions with a propeller angular step of 6 degrees. This results in a total 
simulation time of 1 .0s and 60 time steps per propeller turn or 120 time steps total. 
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In our first simulation, we have selected a mollifier radius Rl D CH =0.01 where D CH =5.1m. Unfortunately, the 
execution of the computer code for this case has automatically been interrupted at t = 0.2 Is , with the error message: 
'a highly deformed element has been detected'. This message has to do with a failure of the adaptive quadrature 
over an element and with the limits (accuracy and maximum number of integration points) set by the user of the 
code. Figure 20 shows the free shear layer geometry at time? = 0.2 Is, which corresponds to the last instance before 
the automatic code interruption. An extensive instability of the shear layer geometry is observed at the junction of 
the starting vortex with the hub vortex, which explains the automatic code interruption with the corresponding 
error message. By increasing the mollifier radius to Rl D CH =0.05, the instability is suppressed allowing thus the 
finishing of the simulation. Figures 21 and 22 show the 3-D and YZ projected views of the free shear layer when 
R / D ch = 0.05 . From the transparent views, it is shown that at the rear downstream part of the free wake some 
instability still occurs. By further increasing the mollifier radius to Rl D CH =0.1, the free wake instability is almost 
fully eliminated in figures 23 and 24. From figures 21 to 24, we observe that, irrespective of the filter used, the shear 
layers always occupy a certain region of the 3-D space with a well-recognized pattern similar to the shape of a 
'mushroom'. Figures 25 and 26 compare slices of the shear layers with planes, which either contain the X axis or are 
normal to it. Using the deformation of the shear layers as a visualizer of the large scale vortices created by the 
propeller, we conclude that the wake of a propeller evolves as a result of a redistribution of both the starting and 
the tip vortices. The starting vortices are mainly responsible for the creation of angular momentum in the propeller 
race, while the tip vortices are connected with the creation of an axial momentum. With the passage of time the 
strong tip vortices, interact with the starting vortices, increasing their diameters and producing a 'mushroom' type 
downstream wake pattern. This mushroom vorticity pattern strengthens with time, accelerating the flow in the 
propeller disk and contracting the flow downstream the propeller disk. Our numerical findings are experimentally 
confirmed in the work of Stettler (2004) [28]. 

Figure 27 presents comparisons of the burst starting propeller thrust forces obtained using either a free WM or a 
frozen WM with experimental results. For the free WM, the following mollifier radii have been used: 

R/D ch =0.01, 0.03, 0.05, 0.07, 0.10 

Curve labeling follows that introduced in the example of a burst starting wing. It is shown that with the exception 
of an initial transient period, the free WM results in practically the same thrust force, irrespective of the selected 
Rl D ch . This force is very near to the experimental measurements, Oosterveld & Oossannen (1975) [15]. On the 
other hand, very large differences exist with the predictions obtained using a frozen WM. At the end of the 
simulation, those differences are of the order of 30% with a trend to increase with the passage of time in figure 27. 

z 

Lx 


i of starting 
tip vortex 


at the junction 
ortex with the 



FIGURE 20. FREE SHEAR LAYER OF A PROPELLER AT t = 0.2 Is . EXTENDED SHEAR LAYER INSTABILITIES OCCUR AT THE JUNCTION 
OF THE HUB VORTEX WITH THE STARTING VORTEX. BLADES COLORED WITH POTENTIAL INTENSITY. SHEAR LAYERS COLORED 

WITH DIPOLE INTENSITY. J = 0.2, R/D CH = 0.01 . 
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FIGURE 25. SLICE OF THE FREE SHEAR LAYER WITH A 
PLANE CONTAINING THE Z AXIS AT t = 1.0s . 
HAIR LINE: R/D CH = 0.05 , BOLD LINE: R/D CH = 0.1 
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FIGURE 27. PROPELLER THRUST DIVIDED BY FLUID DENSITY p AS FUNCTION OF TIME, J = 0.2 . PREDICTIONS WITH DIFFERENT 

WM'S AND COMPARISON WITH EXPERIMENTS. 


Case with J = 0.7 

This case is similar to the previous one in all aspects except velocity of advance which has been increased to 
V = 7. 98m/ s . Again flow simulation lasts for two full propeller revolutions with a propeller angular step of 6 
degrees. This results in a total simulation time of 1 .Os and 60 time steps per propeller turn or 120 time steps total. 

Figure 28 shows the free shear layer geometry for a mollifier radius R! D CH =0.01. Due to the reduced loading 
( J = 0.7 ), instabilities in the free shear layer have been limited at the outmost (downstream) parts of it. By 
increasing the mollifier radius to either R/ D CH =0.05 or Rl D CH =0.10, the instabilities are suppressed in figures 29 
and 30 respectively. From figures 28 to 30, we observe that the 'mushroom' wake pattern remains. Figures 32 and 
33 compare slices of the shear layers with planes which either contain the X axis or are normal to it. 

Figure 31 shows the propeller with a frozen wake model. Comparing figures 28-30 with 31 we observe that, at least 
near the trailing edge, the frozen wake model is a good approximation for the free wake model at this advance 
coefficient. 

Figure 34 presents comparisons of the thrust forces obtained using either a free WM or a frozen WM with 
experimental results. For the free WM, the following mollifier radii have been used: 

R/D ch =0.01, 0.03, 0.05, 0.07, 0.10 

It is shown that the free WM results in practically the same thrust force, irrespective of the selected R / D CH . This 
force is very near to the experimental measurements, Oosterveld & Oossannen (1975) [15]. On the other hand, 
predictions obtained using a frozen WM differ from the experimental measurements by 4% (at the end of the 
simulation) with a trend to increase with the passage of time in figure 34. This explains why the traditional 
generalized wake model used in various propeller theories (lifting line/surface theory and traditional panel 
methods) gives good predictions around the propeller design point. 
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FIGURE 28. 3-D VIEW OF A FREE SHEAR LAYER AT t = l.0s. BLADES COLORED WITH POTENTIAL INTENSITY. SHEAR LAYERS 

COLORED WITH DIPOLE INTENSITY. J = 0.7, R/D CH = 0.01 . 



FIGURE 29. 3-D VIEW OF A FREE SHEAR LAYER AT t = 1 .Os . J = 0.7, R/D CH = 0.05 . 
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FIGURE 30. 3-D VIEW OF A FREE SHEAR LAYER AT t = 1.0s. J = 0.7, R/D CH = 0. 10 . 



FIGURE 31. 3-D VIEW OF A FROZEN WAKE MODEL AT t = 1.0s , / = 0.7 . 
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FIGURE 32. SLICE OF THE FREE SHEAR LAYER WITH A PLANE CONTAINING THE Z AXIS AT t = 1.0s . 
HAIR LINE: R/D CH = 0.01 , BOLD LINE: R/D CH = 0.07 



FIGURE 33. SLICE OF THE FREE SHEAR LAYER AT: t = 1.05, x = -3.36m .HAIR LINE: R/D CH = 0.01 , BOLD LINE: R/D CH = 0.07 



FIGURE 34. PROPELLER THRUST DIVIDED BY FLUID DENSITY p AS FUNCTION OF TIME, J = 0.1 . PREDICTIONS WITH DIFFERENT 

WM'S AND COMPARISON WITH EXPERIMENTS. 
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The Shear Layer Geometry of a Novel Flexible Oscillating Duct Propulsor 

Our last application is the FOD (Flexible Oscillating Duct) propulsor, which conceptually has its origins in the 
natural paradigm of a jellyfish. Jellyfish propulsion mechanics has been investigated by Dabiri (2005, 2006) [8-9], 
who has made many numerical simulations and experimental investigations and has shown the feasibility of using 
variable diameter nozzles as propulsors. In recent years, the development of EAP (Electroactive Polymer) materials 
has opened the way for the construction of bulk artificial muscles, mimicking the motion of aquatic animals, Cohen 
ed. (2004) [5]. A characteristic property for such type of future propulsion systems is their active flexibility through 
an electrically controlled artificial muscle system. The FOD can be considered as a candidate EAP propulsion 
system, since it can produce significant thrust with good hydrodynamic efficiencies. The FOD has been introduced 
for the first time in Politis & Tsarsitalidis (2011) [21], where its feasibility as a ship propulsion mechanism has been 
presented and discussed. 

The FOD is an axisymmetric duct which changes its shape with time in a specific way, determined by the motion of 
a section of it in a radial plane containing the FOD axis. The section of the FOD with this plane has the shape of a 2- 
D airfoil and its motion is similar to that of a biomimetic foil. More specifically, it performs a parallel advancement 
along the FOD axis superimposed to a heaving motion (normal to the axis of the FOD), a pitching motions (with 
angle measured from the axis of the FOD) and an oscillatory change in camber distribution, where all oscillatory 
motions have the same frequency, each with its own phase, Politis & Tsarsitalidis (2011) [21]. 

In the following example, the FOD has a diameter D = 2m , measured at the point of foil section pitch-axis, which is 
situated at 33% chord from the leading edge. The foil section has a chord c = 1.0 m . A NACA 6412 section has been 
used. The FOD advances with velocity V = 4.0 m / s . Each FOD's section performs a heave oscillation with 
amplitude \ =0.5 m, a pitch oscillation with amplitude 25.0 deg and phase 90.0 deg (with respect to heave) and a 
camber oscillation with maximum camber amplitude = m m2 ( m 6412 is the maximum camber of the NACA 6412 
section) and phase -90.0 deg . Oscillation frequency for all motions is 2.0rps . BEM discretization of the FOD has 
been succeeded using 72 elements in the circumferential direction and 16 in the chord-wise direction (same number 
for face and back). Simulation lasts for two periods with 63 time steps per period or 126 time steps total. Total 
simulation time is 1.0s . The Strouhal number for the FOD is St = 0.5 . A pressure type Kutta has been used in the 
simulation. 

Figure 35 shows the free shear layer geometry for a mollifier radius R / D CH = 0.005 where D CH - 2m . Extensive free 
vortex instabilities are observed. By increasing the mollifier radius to Rl D CH =0.10, the instabilities are eliminated 
in figure 36. From the constant dipole lines which are used to color the shear layers in figures 35 and 36, we 
conclude that the shear-layers are composed from a set of circular ring-vortices of variable intensity, coaxial with 
the FOD. Furthermore, irrespective of the filter used, the free shear layers occupy a certain region of the 3-D space 
with a well-recognized pattern. This pattern resembles the shape of a row of 'flying saucers' (one per period) with 
diameters proportional to their age (i.e. older 'flying saucers' have greater diameters). 

Figure 37 presents comparisons of the FOD thrust force obtained using either a free WM or a frozen WM. For the 
free WM, the following mollifier radii have been used: 

R/D ch =0.005, 0.01, 0.05, 0.10 

Regarding predicted thrust force using the free WM, we observe that with the exception of the time neighborhoods 
of t = 0.26s and t = 0.76s , the results are almost independent of the mollifier radius as in all previous examples. For 
the time neighborhoods around t = 0.26s and^ = 0.76s, which coincide with the points where thrust is maximized 
(in absolute value), mollifier independent thrust is obtained only for mollifier radii less than 0.01. Curiously 
enough, around the third points = 0.54s, where thrust is maximized, the calculated thrust is nearly independent 
from the mollifier radius in the whole range of its variation in figure 37. To understand the cause of this behavior 
we concentrate at the structure of the free wake vorticity and more specifically to its relative position with the FOD, 
at points where thrust is maximized. The structure of the free wake vorticity at t = 0.54s and t = 0.7 6s for 
R/ D ch =0.10 is shown in figures 38 and 39 respectively. From figure 38, it is shown that at t = 0.54s the thrust- 
maximum occurs during the FOD expansion phase, while at t = 0.76s the thrust-maximum occurs during the FOD 
contraction phase. Furthermore, the corresponding free wake positions, with respect to the FOD, are entirely 
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different for the two cases; it is a fact which explains why the value of the mollifier radius affects differently the 
calculated thrust. 

More specifically, during the expansion phase, the stronger wake ring vortices are relatively far from the FOD 
section. As a result, the calculated thrust is independent from the mollifier radius in the whole range of its variation. 
During the FOD contraction phase, the stronger wake ring- vortices are closer to the FOD section trailing edge. This 
has a pronounced effect on the Kutta condition and consequently to the calculated FOD pressures and forces. The 
effect of the mollifier radius to the position of the free shear layer at t = 0.16s is shown in figure 40, where hair lines 
denote the free shear layer geometry for R / D CH = 0.005 while bold lines denote free shear layer geometry for 
R/ D ch =0.10. From this figure, we observe that the mollifier radius affects the position of the shear layer in the 
vicinity of the FOD trailing edge, explaining thus the dependence of thrust from the mollifier radii at t = 0.16s as 
shown in figure 37. 

Finally notice that in the case of a FOD, thrust predictions obtained using a frozen WM differ considerably from 
that with a free WM in figure 37. 



FIGURE 35. FREE SHEAR LAYER OF A FOD AT t = 1.0s. EXTENSIVE SHEAR LAYER INSTABILITIES OCCUR. FOD COLORED WITH 
POTENTIAL INTENSITY. SHEAR LAYER COLORED WITH DIPOLE INTENSITY. St = 0.5, R/D CH = 0.005 



FIGURE 36. FREE SHEAR LAYER OF A FOD AT t = 1 .0s . VORTEX INSTABILITIES HAVE BEEN SUPPRESSED. St = 0.5, R / D CH =0.10 
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FIGURE 37. THRUST (X-AXIS) FORCE OF THE FOD - PREDICTIONS WITH DIFFERENT WAKE MODELS. 



FIGURE 38. FREE SHEAR LAYER OF A FOD AT t = 0.54s . R/D CH = 0.10 



FIGURE 39. FREE SHEAR LAYER OF A FOD AT t = 0.76s . R/D CH = 0.10 


25 


www.daoe-journal.org 


Development and Applications of Oceanic Engineering (DAOE) Volume 5 2016 



FIGURE 40. SLICE OF THE FREE SHEAR LAYER DURING FOD CONTRACTION ( t = 0.16s, z = 0.0 m ) 
HAIR LINE: R/D CH = 0.005 , BOLD LINE: R/D CH = 0. 10 


Conclusions 

The details of a filtering technique are presented, which in combination with a boundary element formulation and 
a time stepping algorithm can be used to simulate the motion of the free shear layers, produced in cases of 
incompressible non-viscous flows around systems of unsteadily moving rigid or flexible lifting bodies. This is a 
major innovation in treating complex unsteady propulsion problems since no simplifying assumptions (like that of 
a frozen or generalized wake) are necessary. The method has been applied to a number of flows around lift 
producing configurations, performing unsteady motions with practical interest. From the systematic calculations, 
the following conclusions can be drawn: 

In all cases, there is a certain non-dimensional range for the mollifier radius, almost problem independent, where 
the free wake vorticity is trapped inside a subset of the 3-D space and forms a well-recognized pattern. For the 
simulated cases, this 'attracting configuration' varies considerably from case to case. More specifically: (a) for a 
steadily moving wing it is a quadrilateral ring-vortex with roll-up at its free edges (tip vortices and starting vortex), 
(b) in the case of a biomimetic wing it is a repeated S-shaped pattern of oblique ring vortices or equivalently a 
reverse Karman Vortex Street, (c) in the case of a burst starting propeller it is a 'mushroom' type 3-D vortex pattern, 
while (d) in the case of a FOD it has the shape of a row of 'flying saucers' with increasing diameters, the number of 
which is equal to the number of the simulation periods. 

For smaller values of the mollifier, radius instabilities or even chaotic behavior are almost always present in the 
distribution of the free vorticity inside the attracting configuration. By increasing the mollifier radii or using age 
dependent mollifier radii, those instabilities gradually disappear. It is of great importance that chaotic distributions 
of the wake vorticity, obtained for the lower mollifier radii, result in predictions for forces very close to each other, 
as far as the vorticity remains inside its attracting configuration and it is not so close to the body boundaries. This 
result indicates that under the seemingly chaotic picture in the free shear layer vortices, some type of order still 
exists. 

From the simulations, it is shown that in the range of the Strouhal numbers where natural flyers and swimmers 
operate, the frozen wake model can be an attractive alternative for preliminary force predictions. 
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